Minimum Spanning Networks

Purpose
Required packages and data
Full MSN
MSN per year
ZONE 1
ZONE 2
By ZONE2 and Year

Purpose

This document will explore the distribution of MLGs across years as well as filtering strategies.

Required packages and data

library(PramCurry)

## Loading required package: poppr
## Loading required package: adegenet
## Loading required package: ade4
##    ==========================
##     adegenet 1.4-2 is loaded
##    ==========================
## 
##  - to start, type '?adegenet'
##  - to browse adegenet website, type 'adegenetWeb()'
##  - to post questions/comments: adegenet-forum@lists.r-forge.r-project.org
## 
## 
## This is poppr version 1.1.2.99.70. To get started, type package?poppr

library(reshape2)
library(ggplot2)
library(ape)
library(poppr)
library(adegenet)
library(igraph)

## 
## Attaching package: 'igraph'
## 
## The following object is masked from 'package:ape':
## 
##     edges

library(animation)
options(stringsAsFactors = FALSE)
data(ramdat)
data(pop_data)
data(myPal)
sessionInfo()

## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] animation_2.3     igraph_0.7.1      ape_3.1-4.5       ggplot2_1.0.0    
## [5] reshape2_1.4      PramCurry_1.0     poppr_1.1.2.99-70 adegenet_1.4-2   
## [9] ade4_1.6-2       
## 
## loaded via a namespace (and not attached):
##  [1] assertthat_0.1      bitops_1.0-6        boot_1.3-11        
##  [4] caTools_1.17.1      colorspace_1.2-4    digest_0.6.4       
##  [7] dplyr_0.2           evaluate_0.5.5      fastmatch_1.0-4    
## [10] formatR_1.0         gdata_2.13.3        ggmap_2.3          
## [13] grid_3.1.2          gtable_0.1.2        gtools_3.4.1       
## [16] htmltools_0.2.6     httpuv_1.3.0        knitr_1.6          
## [19] labtools_0.4.1      lattice_0.20-29     lubridate_1.3.3    
## [22] mapproj_1.2-2       maps_2.3-9          MASS_7.3-34        
## [25] Matrix_1.1-4        memoise_0.2.1       munsell_0.4.2      
## [28] nlme_3.1-117        parallel_3.1.2      pegas_0.6          
## [31] permute_0.8-3       phangorn_1.99-7     plotrix_3.5-7      
## [34] plyr_1.8.1          png_0.1-7           proto_0.3-10       
## [37] Rcpp_0.11.2         RgoogleMaps_1.2.0.6 rjson_0.2.14       
## [40] RJSONIO_1.3-0       rmarkdown_0.3.10    scales_0.2.4       
## [43] shiny_0.10.1        stringr_0.6.2       tools_3.1.2        
## [46] vegan_2.0-10        xtable_1.7-4        yaml_2.1.13

Full MSN

First thing to do is to create the full MSN across years. Note that the get_layout function will preserve the layout (“MASTER”) so that it can be subsetted.

setpop(ramdat) <- ~ZONE2
options(digits = 10)
(newReps <- other(ramdat)$REPLEN)

##  PrMS6A1  Pr9C3A1 PrMS39A1 PrMS45A1 PrMS43A1 
##        3        2        2        4        4

newReps[3] <- 4
(newReps <- fix_replen(ramdat, newReps))

##  PrMS6A1  Pr9C3A1 PrMS39A1 PrMS45A1 PrMS43A1 
##  3.00000  2.00000  4.00001  4.00000  4.00000

b.msn.full <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
                        include.ties = TRUE)

goodSeed <- 19
goodSeed2 <- 17
goodSeedAlt <- 14
# for (i in goodSeed:100){
  set.seed(goodSeed)
#   set.seed(i)
  MASTER <- get_layout(b.msn.full$graph)
  plot_poppr_msn(ramdat, b.msn.full, gad = 10, palette = funky, mlg = TRUE, 
                 layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
                 quantiles = FALSE,
                 vertex.label.color = "firebrick", inds = "none")

plot of chunk unnamed-chunk-2

#   prompt <- paste("save", i, "as seed?")
#   accept <- readline(prompt)
#   if (substr(accept, 1, 1) == "y"){
#     theSeed <- i
#     break
#   }
# }

MSN per year

Now we will see how the MSNs are plotted per year. For this, we will do two things:

color the nodes based on the full MSN
ensure that the node coordinates are equal

saveGIF({
  plot_poppr_msn(ramdat, b.msn.full, gad = 10, palette = funky, mlg = TRUE, 
                   layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
                   quantiles = FALSE,
                   vertex.label.color = "firebrick")
  
  for (i in levels(pop(ramdat))){
    population_subgraph(ramdat, b.msn.full, 
                        quantiles = FALSE,
                        gad = 10, 
                        palette = funky, 
                        mlg = TRUE, 
                        layfun = MASTER, 
                        nodebase = 1.75, 
                        vertex.label.font = 2,
                        cutoff = 0.148,
                        nodelab = 1000,
                        vertex.label.color = "firebrick",
                        sublist = i)
  }
}, movie.name = "by_year.gif", interval = 1, ani.height = 1000, ani.width = 1000)

## Executing: 
## 'convert' -loop 0 -delay 100 Rplot1.png Rplot2.png Rplot3.png
##     Rplot4.png Rplot5.png Rplot6.png Rplot7.png Rplot8.png
##     'by_year.gif'
## Output at: by_year.gif

## [1] FALSE

# pdf("by_year.pdf", width = 10, height = 10)
plot_poppr_msn(ramdat, b.msn.full, gad = 10, palette = funky, mlg = TRUE, 
                 layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
                 quantiles = FALSE,
                 vertex.label.color = "firebrick")

plot of chunk unnamed-chunk-3

for (i in levels(pop(ramdat))){
  population_subgraph(ramdat, b.msn.full, 
                      quantiles = FALSE,
                      gad = 10, 
                      palette = funky, 
                      mlg = TRUE, 
                      layfun = MASTER, 
                      nodebase = 1.75, 
                      vertex.label.font = 2,
                      cutoff = 0.148,
                      nodelab = 1000,
                      vertex.label.color = "firebrick",
                      sublist = i)
}

plot of chunk unnamed-chunk-3

# dev.off()

# Just to see once more with only allowing 1 mutational step for Bruvo's distance
plot_poppr_msn(ramdat, b.msn.full,
               quantiles = FALSE,
               gad = 10, 
               palette = funky,
               mlg = TRUE, 
               layfun = MASTER, 
               nodebase = 1.75, 
               vertex.label.font = 2,
               cutoff = 0.051,
               vertex.label.color = "firebrick")

plot of chunk unnamed-chunk-3

ZONE 1

setpop(ramdat) <- ~ZONE1
# zonePal <- get_year_pal(ramdat)
b.msn.zone1 <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
                        include.ties = TRUE)

plot_poppr_msn(ramdat, 
               b.msn.zone1, 
               gadj = 10, 
               palette = funky, 
               mlg = TRUE, 
               layfun = MASTER, 
               nodebase = 1.75, 
               vertex.label.font = 2,
               quantiles = FALSE,
               vertex.label.color = "firebrick")

plot of chunk unnamed-chunk-5

for (i in levels(pop(ramdat))){
  population_subgraph(ramdat, 
                      b.msn.zone1, 
                      quantiles = FALSE,
                      gad = 10, 
                      palette = funky, 
                      mlg = TRUE, 
                      layfun = MASTER,
                      nodebase = 1.75, 
                      vertex.label.font = 2,
                      cutoff = 0.148,
                      vertex.label.color = "firebrick",
                      sublist = i,
                      nodelab = 1000,
                      inds = "none")
}

plot of chunk unnamed-chunk-5

ZONE 2

setpop(ramdat) <- ~ZONE2
# zonePal <- get_year_pal(ramdat)
b.msn.zone <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
                        include.ties = TRUE)

zdegs <- igraph::degree(b.msn.zone$graph)
names(zdegs) <- igraph::V(b.msn.zone$graph)$label
sort(zdegs)

## MLG.26 MLG.19 MLG.65 MLG.64 MLG.55 MLG.58 MLG.43  MLG.7 MLG.40 MLG.11 
##      1      1      1      1      1      1      1      1      1      1 
## MLG.42 MLG.38  MLG.9 MLG.62 MLG.10 MLG.37 MLG.31 MLG.61  MLG.1 MLG.63 
##      1      1      1      1      1      1      1      1      1      1 
## MLG.20 MLG.69 MLG.67 MLG.70 MLG.32 MLG.25  MLG.5 MLG.44 MLG.36 MLG.35 
##      2      2      2      2      2      2      2      2      2      2 
## MLG.13  MLG.2 MLG.45  MLG.8 MLG.27 MLG.16 MLG.41  MLG.3 MLG.51 MLG.66 
##      2      2      2      2      3      3      3      3      3      3 
## MLG.33 MLG.59  MLG.6 MLG.57 MLG.56 MLG.50 MLG.34 MLG.39  MLG.4 MLG.24 
##      3      3      3      3      3      3      3      3      4      4 
## MLG.28 MLG.68 MLG.60 MLG.52 MLG.46 MLG.54 MLG.15 MLG.47 MLG.21 MLG.29 
##      4      4      4      4      4      4      4      4      5      5 
## MLG.14 MLG.30 MLG.18 MLG.12 MLG.49 MLG.23 MLG.53 MLG.48 MLG.22 MLG.17 
##      5      5      5      5      5      6      6      6      7      8

edges_to_remove <- E(b.msn.zone$graph)[E(b.msn.zone$graph)$weight > 0.05]
clusts <- clusters(delete.edges(b.msn.zone$graph, edges_to_remove))
(names(clusts$membership) <- V(b.msn.zone$graph)$label)

##  [1] "MLG.22" "MLG.27" "MLG.26" "MLG.4"  "MLG.24" "MLG.21" "MLG.16"
##  [8] "MLG.20" "MLG.19" "MLG.41" "MLG.17" "MLG.28" "MLG.29" "MLG.23"
## [15] "MLG.68" "MLG.69" "MLG.3"  "MLG.51" "MLG.14" "MLG.65" "MLG.64"
## [22] "MLG.67" "MLG.70" "MLG.60" "MLG.52" "MLG.66" "MLG.30" "MLG.46"
## [29] "MLG.33" "MLG.32" "MLG.59" "MLG.55" "MLG.58" "MLG.53" "MLG.18"
## [36] "MLG.54" "MLG.43" "MLG.12" "MLG.48" "MLG.7"  "MLG.25" "MLG.5" 
## [43] "MLG.40" "MLG.11" "MLG.6"  "MLG.44" "MLG.42" "MLG.38" "MLG.9" 
## [50] "MLG.62" "MLG.36" "MLG.15" "MLG.10" "MLG.37" "MLG.35" "MLG.31"
## [57] "MLG.13" "MLG.57" "MLG.61" "MLG.47" "MLG.1"  "MLG.49" "MLG.56"
## [64] "MLG.2"  "MLG.50" "MLG.63" "MLG.34" "MLG.45" "MLG.39" "MLG.8"

clustPal  <- clusts$csize
theClusts <- clustPal > 3
clustOrder <- order(clustPal, decreasing = TRUE)
clustPal[clustOrder][theClusts[clustOrder]]   <- RColorBrewer::brewer.pal(sum(theClusts), "Set1")
clustPal[clustOrder][!theClusts[clustOrder]]  <- gray.colors(sum(!theClusts))
nodeList <- lapply(1:length(clustPal), function(x) which(clusts$membership == x))


# saveGIF({
  
  plot_poppr_msn(ramdat, b.msn.zone, gad = 10, palette = funky, mlg = TRUE, 
                 layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
                 quantiles = FALSE,
                 vertex.label.color = "firebrick",
                 mark.groups = nodeList[theClusts], 
                 mark.border = clustPal[theClusts],
                 mark.col = transp(clustPal[theClusts], 0.05),
                 mark.expand = 2,
                 mark.shape = 0)

$plot of chunk unnamed-chunk-6$

  for (i in levels(pop(ramdat))){
    population_subgraph(ramdat, 
                        b.msn.zone, 
                        gadj = 10,
                        palette = funky, 
                        mlg = TRUE,
                        layfun = MASTER,
                        nodebase = 1.75,
                        vertex.label.font = 2,
                        quantiles = FALSE,
                        vertex.label.color = "firebrick",
                        sublist = i, 
                        inds = "none",
                        nodelab = 1000)
  }

plot of chunk unnamed-chunk-6

# }, movie.name = "by_zone2.gif", interval = 1, ani.height = 1000, ani.width = 1000)

mainMLGs <- order(table(ramdat@mlg), decreasing = TRUE)[1:10]

pdf("by_zone2.pdf", width = 352/2.54/10, height = 352/2.54/10, pointsize = 40)
plot_poppr_msn(ramdat, b.msn.zone, gad = 10, palette = funky, mlg = TRUE, 
               layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
               quantiles = FALSE, 
#                inds = mainMLGs, nodelab = 1000,
               vertex.label.color = "firebrick",
               mark.groups = nodeList[theClusts], 
               mark.border = clustPal[theClusts],
               mark.col = transp(clustPal[theClusts], 0.05),
               mark.expand = 2,
               mark.shape = 0)
for (i in levels(pop(ramdat))){
  population_subgraph(ramdat, 
                      b.msn.zone, 
                      gadj = 10,
                      palette = funky, 
                      mlg = TRUE,
                      layfun = MASTER,
                      nodebase = 1.75,
                      vertex.label.font = 2,
                      quantiles = FALSE,
                      vertex.label.color = "firebrick",
                      sublist = i, 
                      inds = "none",
                      nodelab = 1000)
}
dev.off()

## pdf 
##   2

By ZONE2 and Year

zone2pal <- char2pal(setpop(ramdat, ~ZONE2)@pop.names, funky)
setpop(ramdat) <- ~ZONE2/Pop
zone2yearpal <- zone2pal[sub("_.+?$", "", ramdat@pop.names)]
names(zone2yearpal) <- ramdat@pop.names
# zonePal <- get_year_pal(ramdat)
b.msn.zone.year <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
                        include.ties = TRUE)
# saveGIF({
# pdf("by_zone2year.pdf", width = 10, height = 10)
  plot_poppr_msn(ramdat, b.msn.zone.year, gad = 10, palette = funky, mlg = TRUE, 
               layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
               quantiles = FALSE,
               vertex.label.color = "firebrick")

plot of chunk unnamed-chunk-7

  for (i in sort(levels(pop(ramdat)))){
    population_subgraph(ramdat, 
                        b.msn.zone.year, 
                        gadj = 10,
                        mlg = TRUE,
                        layfun = MASTER,
                        nodebase = 1.75,
                        vertex.label.font = 2,
                        quantiles = FALSE,
                        palette = funky,
                        vertex.label.color = "firebrick",
                        sublist = i, 
                        inds = "none",
                        nodelab = 1000,
                        main = paste("\n\n", i))
  }

plot of chunk unnamed-chunk-7

# dev.off()
# }, movie.name = "by_zone2year.gif", 
# interval = 1, ani.height = 1000, ani.width = 1000)

Minimum Spanning Networks

Zhian N. Kamvar

September 29, 2014

Purpose

Required packages and data

Full MSN

MSN per year

ZONE 1

ZONE 2

By ZONE2 and Year